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1.  INTRODUCTION 


Most  processes  occurring  on  a  crystal  surface  and  at  interfaces  are  expected  to  be  influenced  by  the 
kinetics  of  the  atoms.  The  mobility  of  atoms  on  a  crystal  surface  can  be  orders  of  magnitude  larger  than 
in  the  bulk  at  the  same  temperature  and  become  the  determining  step  in  many  processes  such  as  surface 
reconstruction,  chemical  reaction,  and  epitaxial  growth.  Diffusion  on  metals  continues  to  receive  both 
experimental  and  theoretical  attention  (Ehrlich  and  Stolt  1980;  Doll  and  Voter  1987;  Bonzel  1983;  Ehrlich 
1990).  The  relation  between  the  surface  diffusivity  and  surface  structure  has  been  studied  previously  by 
detailed  molecular  dynamic  (MD)  simulations  using  simple  model  potential  energy  functions,  which  are 
typically  pairwise  additive.  These  studies  (De  Lorenzi  1986;  De  Lorenzi  and  Jacucci  1985)  have  shown 
that  the  surface  structure  and  the  softness  of  the  repulsive  potential  play  a  dominant  role  in  determining 
the  diffusion  properties  and  that  the  migration  of  adatoms  can  occur  not  only  through  simple  hops  between 
binding  sites  on  a  surface,  but  through  more  complicated  mechanisms,  such  as  adatom  exchange  on 
channelled  surfaces.  In  modelling  metallic  materials,  however,  simple  pair  potentials  do  not  adequately 
describe  the  elastic  properties  of  the  metal  and  require  a  volume-dependent  energy  term  to  correct  for  this 
inadequacy  (Lee  1981;  Johnson  1988).  Ambiguities  arise  in  defining  the  termination  of  the  volume  at 
surfaces  and  bulk  defects,  thus  introducing  additional  uncertainty  into  the  model. 

The  embedded  atom  method  (EAM)  was  recently  developed  (Daw  and  Baskes  1983, 1984;  Finnis  and 
Sinclair  1981)  as  a  means  of  calculating  realistic  ground-state  properties  of  metallic  systems,  and 
circumvents  the  definition  of  volume  by  using  an  energy  term  which  is  electron-density-dependent,  a 
property  which  is  always  definable.  The  electron-density-dependent  term  is  the  embedding  energy,  or  the 
energy  required  to  embed  an  atom  into  the  local  electron  density  of  all  surrounding  atoms.  To  our 
knowledge,  our  analysis  of  various  embedded-atom-type  modeling  schemes  for  adatom  diffusion  studies 
represents  the  first  study  based  on  the  realistic  bulk  models.  The  major  goal  of  this  study  is  to  examine 
the  use  of  commonly  used  EAM  potentials,  which  have  been  parameterized  to  reproduce  bulk  properties 
in  adatom/surface  studies,  and  to  determine  the  applicability  of  these  potential  energy  functions  for  studies 
of  interfacial  processes.  As  an  initial  step  towards  such  studies,  self-diffusion  on  nickel  surfaces  is  chosen 
as  a  prototype,  but  the  results  can  be  extended  to  any  face-centered  cubic  (FCC)  system. 

For  estimating  the  diffusion  coefficients  and  Arrhenius  parameters,  we  used  the  variational  transition 
state  theory  (VTST)  approach,  instead  of  adapting  the  cost-intensive  MD  simulations,  which  yield  quite 
detailed  information.  VTST  has  been  used  successfully  in  calculating  rates  of  gas  phase  reactions. 
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Additionally,  it  has  been  applied  to  other  gas-surface  diffusion  studies.  We  refrained  from  employing 
molecular  statics  due  to  both  the  anharmonic  features  of  the  potential  energy  surface  and  the  complicated 
saddle  point  surface  in  configuration  space.  In  the  present  study,  we  present  results  of  VTST  calculations 
of  diffusion  of  atomic  nickel  on  several  nickel  surfaces  using  four  different  EAM  potential  functions. 
These  pair  functionals  provide  similar  descriptions  of  the  atomic  interaction,  but  have  completely  different 
functional  forms  for  the  pairwise  interaction  as  well  as  for  the  embedding  function  and  the  atomic  density. 

The  remainder  of  this  report  is  organized  as  follows:  Section  2  provides  details  of  (a)  a  brief  review 
of  the  dynamical  method,  VTST;  (b)  the  background  and  a  short  description  of  the  EAM  potentials  used 
in  this  study;  and  (c)  the  description  of  the  structural  models  of  the  nickel  lattices  used  in  the  calculations. 
Discussion  of  the  results  of  surface  self-diffusion  and  the  properties  of  the  nickel  adatom  on  different 
nickel  surfaces  are  presented  in  Section  3  and  the  results  are  summarized  in  Section  4. 

2.  THEORY  AND  MODELS 

2.1  VTST  Theory.  Surface  diffusion  is  typically  thought  of  as  a  random  walk  model  of  uncorrelated 
hops  of  an  isolated  adatom  to  equivalent  and  adjacent  binding  sites.  Each  hop  is  assumed  to  be  an 
individual  unimolecular  process  with  thermal  rate  constant,  k(T).  The  thermal  diffusion  coefficient,  D(T), 
can  be  described  in  terms  of  the  unimolecular  rate  constant: 

D(T)  -  ilk(T)  (1) 

2y 

where  a  is  the  hop  length  and  y  is  the  dimensionality  of  the  system  (y  =  2  for  diffusion  on  a  surface 
wherein  all  binding  sites  are  equivalent;  and  y  =  1  for  surface  diffusion  along  channels  [Voter  and  Doll 
1984]  such  as  along  the  [110]  or  [001]  directions  on  Ni[l  TO]). 

The  rate  constant  k(T)  is  calculated  from  canonical  VTST  with  semiclassical  adiabatic  ground-state 
transmission  coefficients.  VTST  has  been  used  successfully  in  calculating  such  rates  in  other  surface 
diffusion  studies;  details  are  given  elsewhere  (Truhlar,  Isaacson,  and  Garrett  1983;  Lauderdale  and  Truhlar 
1985;  Truong  and  Truhlar  1987,  1988).  In  review,  the  procedure  consists  of  first  locating  the  reactant, 
product,  and  saddle  point  geometries  (locations  at  which  the  gradient  vanishes),  then  determining  the 
reaction  path  b,  following  the  negative  of  the  gradient  vector  in  mass  weighted  coordinates  from  the 
saddle  point  down  into  the  reactant  and  product  geometries.  The  reaction  coordinate  is  defined  as  the 
distance  on  the  reaction  path  from  the  saddle  point  and,  by  convention,  s  is  negative  on  the  reactant  side 
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of  the  saddle  point.  The  generalized  transition  state  theory  rate  is  expressed  in  terms  of  the  location  of 
the  transition  state  along  the  reaction  coordinate 


where  kB  is  Boltzmann’s  constant;  T  is  the  temperature  of  the  system;  QGT  (T,  s)  is  the  partition  function 
for  the  bound  degrees  of  freedom  at  the  generalized  transition  state  location  at  s;  VMEP  (s)  is  the  potential 
along  the  reaction  path;  h  is  Planck’s  constant;  and  QR  (T)  is  the  partition  function  of  the  reactant  species. 
The  partition  functions  are  defined  with  their  zeros  of  energies  at  the  local  minimum  of  the  potential  for 
the  bound  degrees  of  freedom,  and  can  be  evaluated  from  the  potential  energy  surface  using  standard 
statistical  mechanical  methods.  In  canonical  variational  transition  state  theory,  the  location  of  s  is 
optimized  to  minimize  the  rate.  Quantum  mechanical  effects  such  as  tunneling  ana  zero-point  energy 
effects  are  treated  in  the  VTST  calculation;  however,  these  effects  are  not  important  in  the  systems  studied. 
Additionally,  variationally  optimizing  the  location  of  the  dividing  surface  (i.e.,  the  choice  of  s  in  Equation 
2)  always  yielded  the  saddle  point  (i.e.,  s  =  0). 


2.2  Embedded-Atom  Pair  Functionals.  The  derivation  of  the  embedded-atom  method  and  its 
applications  to  metallic  systems  are  given  elsewhere  (Daw  and  Baskes  1983,  1984;  Finnis  and  Sinclair 
1981).  Its  basis  on  density  functional  theory,  low  computational  demands,  and  success  in  describing  bulk 
phenomena  make  the  method  a  desirable  one  for  use  in  computer  simulation  of  interfacial  processes  in 
metals. 


In  the  EAM,  the  energy  of  a  collection  of  atoms  is  expressed  as  a  sum  of  a  density-dependent 
embedding  energy  and  a  pair  interaction  term, 


u-EwlEW. 

i  L  j*i 

where  Fj  (p()  is  the  energy  needed  to  embed  atom  i  into  the  local  electron  density  p;  of  the  surrounding 
atoms,  and  is  a  pair  potential  between  atoms  i  and  j.  The  functional  foim  of  the  EAM  has  been 
derived  by  Manninen  (1986)  and  Jacobsen,  Nprskov,  and  Puska  (1987)  using  density  functional  theory. 
The  two-body  potential  within  this  effective  medium  theory  arises  from  atomic-core-electron  gas 
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electrostatic  attraction  whereas  it  is  considered  to  be  repulsive  in  some  of  the  empirical  procedures 
described  below.  Further,  in  some  empirical  approaches,  the  two-body  potential  is  separated  into  a  product 
of  two  one-atom  functionals,  analogous  to  the  embedding  energy  term,  while  some  other  EAM-type 
models  have  used  nonseparable  forms.  Different  empirical  schemes  have  thus  been  used  to  construct 
embedded  atom  pair  functionals  for  nickel;  they  are  summarized  below. 

In  die  Foiles,  Baskes,  and  Daw  (FBD)  (1986)  approach,  the  local  electron  density  into  which  atom 
i  is  embedded  is  approximated  by  the  sum  over  atomic  electron  densities  p*  (r-)  from  all  other  atoms  in 
the  system: 

Pi  *  E  Pj  <rij>  •  (4) 

j 

This  functional  form  has  the  advantage  of  including  many-body  interactions  due  to  the  nonlinear  nature 
of  the  embedding  function  even  though  the  computational  effort  of  evaluating  the  embedding  energy  is 
similar  to  that  for  the  pairwise-additive  interactions.  The  functional  form  of  the  pairwise  additive  term 
is: 


OijCr)  «  [Zj(r)Zj(r)  -  Z^Z^r c)]/r  .  r<rc 

-0  r  £  rc  (5) 

where  2^(r)  is  the  effective  charge  of  atom  m,  and  re  is  chosen  to  be  large  enough  that  the  pair 
interactions  are  very  small  at  that  value.  The  effective  charge  of  an  atom  is  given  the  functional  form  as 
follows: 

Z(r)  «  ZqU  +  Prv)exp(-ar)  .  (6) 

The  atomic  densities  are  approximated  by  the  single-determinant  Hartree-Fock  calculations  of  Dementi 
and  Roetti  (1974).  For  nickel,  the  atomic  density  is  further  approximated  as  the  sum  of  electron  densities 
from  the  3d  and  4s  shells,  weighted  by  the  number  of  electrons  in  each  shell.  The  number  of  4s  electrons, 
n,,  is  used  as  an  adjustable  parameter  and  is  set  to  two  for  all  the  calculations  reported  here. 

For  a  given  set  of  parameters  for  the  atomic  densities  and  pair  interaction  terms,  the  embedding 
function  was  then  adjusted  (Foiles  198S)  so  that  the  simplified  equation  of  state  of  the  metal  (Rose,  Smith, 
Guinea,  and  Ferrante  1984)  was  exactly  reproduced  for  successive  values  of  the  density.  The  embedding 
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function  for  nickel  was  represented  in  a  tabular  form  as  a  function  of  atomic  density.  The  four  parameters 
appearing  in  Equation  6  and  the  embedding  function  were  fitted  (Foiles,  Baskes,  and  Daw  1986)  to 
experimental  data  for  the  bulk  (equation  of  state,  elastic  modulii,  defect  energetics,  and  heats  of  solutions 
of  binary  alloys)  by  a  least-squares  procedure  and  those  defining  the  pair  interactions  for  nickel  are  given 
in  Table  1.  In  the  present  application,  the  tabulated  embedding  function  is  fitted  to  a  cubic  spline. 

The  differences  between  the  FBD  and  other  EAM  approaches  are  primarily  found  in  the  derivation 
and  interpretation  of  the  electron  density  surrounding  the  atom  and  the  embedding  function.  Voter  and 
Chen  (VC)  (1987)  have  basically  followed  the  FBD  scheme  described  in  Equations  3  and  4,  with  the 
following  deviations:  The  pairwise  additive  interaction  is  taken  to  be  a  Morse  potential, 

^ij(r)  *  Dm  { 1  -  exp[  -aM(r  -  RM)]}2  -  DM  ,  (7) 

where  Rm,  Dm,  and  aM  define  the  location,  depth,  and  curvature  of  the  minimum,  respectively.  The 
density  function  p(r)  is 

p(r)  *  r6[exp( -(Jr)  +  29exp(-2Pr)]  ,  (8) 

which  is  a  modification  of  the  density  of  a  hydrogenic  4s  orbital.  The  functions  d>(r)  and  p(r)  and  the 
corresponding  first  derivatives  are  forced  to  go  smoothly  to  zero  at  a  cutoff  distance,  rc.  The  four 
parameters  for  Equations  7  and  8  at  the  cutoff  rc  and  the  embedding  function  were  fitted  to  the  same 
experimental  data  used  in  the  FBD  fit  (Table  1).  Like  the  FBD  embedding  function,  the  VC  embedding 
function  is  obtained  in  terms  of  tabular  form. 

Oh  and  Johnson’s  (OJ)  (1988)  approach  was  to  develop  an  EAM  model  with  simple  functional  forms 
which  result  in  fewer  parameters  to  fit  and  thus  provide  easy  access  for  adapting  the  EAM  model.  The 
embedding  energy  F(p)  is  defined  as 

F(p)  -  o(Pi/Pe)n  ♦  b(Pi/Pc)  ,  (9) 

where 

Pi<r)  *  E  f<r>  •  00) 
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Table  1.  Parameters  for  the  EAM  Potentials 


FBD  Potential 

VC  Potential 

Zo(eV-Ar1/2 

37.9326 

DM(eV) 

1.5335 

a  (A"1) 

1.8633 

«m  (A-1) 

1.7728 

P  (A-v) 

0.8957 

Pm  (A'1) 

2.2053 

V 

1 

Rm  (A) 

3.6408 

NS 

2 

rcut  (A) 

4.7895 

rcut  (A) 

4.8 

|  OJ  Potential* 

ATVF  Potential11 

P 

6 

ai 

29.057085 

Y 

9 

h. 

-76.04625 

5 

20 

23 

48.08920 

rcut 

1.9 

a4 

-25.96604 

fe 

1.0 

a5 

79.15121 

♦e 

0.63199 

rl 

1.2247449 

a 

-3.1684 

h 

1.1547054 

b 

-5.1237 

r3 

1.1180065 

n 

0.12501 

r4 

1.0000000 

Pe 

12.551 

r5 

0.8660254 

— 

— 

A1 

60.537985 

— 

— 

a2 

-80.102414 

— 

— 

Ri 

1.2247449 

— 

R2 

1.1180065 

a  JJ4 

0e  and  r^,,  are  expressed  in  eV  and  nearest  neighbor  distance  *^="  A,  respectively. 
The  other  parameters  are  dimensionless. 

“The  coefficients  are  in  units  of  eV  and  the  lattice  parameter,  3.52  A,  respectively. 
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and 

f(r)  -  feexp(-p[r/re  -  1  ])  -  feexp(  -  p  [rc/re  -  1  ]) 

-  P f c (  1  "  exp(8[r/rc  -  rc/re])}  exp( -p[rc/re  -  1  ])/5  .  (11) 


The  pairwise  additive  interaction  has  a  similar  form 

O(r)  «  Oe exp ( —  Y  [ r / re  -  1])  -  <Dcexp(-Y  [rc/re  -  1]) 

-y<De  {1  -  exp  (8  [r/re-rc/re]))  exp  (~y  trc/re  -  1  ])/8  .  (12) 


The  nine  parameters  for  both  embedding  and  pair  potential  terms  (Equations  9-12)  arc  given  in  Table  1. 

Ackland,  Tichy,  Vitek,  and  Finnis  (ATVF)  (1987)  identified  F(p)  with  the  second  moment  of  the 
density  of  states  and  thus  used  a  square-root  function  for  F(p).  The  function  for  both  F(p)  and  <X>  consist 
of  a  series  expansion  with  built-in  step  functions  (denoted  H(x)  in  Equation  17),  the  coefficients  of  which 
have  been  fitted  to  the  experimental  data  as  well  as  the  pressure-volume  values  calculated  from  first- 
principles.  The  forms  of  the  functions  are  as  follows: 


F(Pi)  -  -  ,  (13) 

where 

Pi  -  E  ♦  (rtj)  04) 

and 

2 

4»  (r)  -  £  (Ak(Rk  -  r )3  H  (Rk  -  r)J  .  (15) 

k-1 


The  pairwise  additive  term  <t>  has  the  form 

5 

*<r)»  £  {ak(rk-r)3H(rk-r)}  .  (16) 

k*  1 
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The  step  function  H(x)  used  in  Equations  IS  and  16  has  the  form 


H(x)  =  0,  x  <  0 

=  1,  x  £  0.  (17) 

The  14  parameters  which  include  the  embedding  energy  term  for  the  ATVF  potential  energy  function  are 
given  in  Table  1. 

2.3  Lattice  Models.  The  structural  model  used  in  self-diffusion  on  Ni(100)  consists  of  162  atoms 
over  4  layers;  41  atoms  are  in  each  of  the  first  and  third  layers,  and  40  atoms  are  in  each  of  the  second 
and  fourth  layers.  Figure  la  illustrates  the  crystal  model  used  in  this  calculation;  the  six  shaded  atoms 
denote  the  first-layer  atoms  which  were  allowed  to  move  in  this  calculation.  The  reactant  and  product 
sites  are  labelled  R  and  P,  respectively,  and  lie  over  the  two  second-layer  atoms  which  were  allowed  to 
move  in  this  calculation.  The  transition  state,  which  is  located  directly  between  the  reactant  and  product 
sites,  is  labelled  T. 

The  model  of  the  Ni(l  1 1)  crystal  consists  of  120  nickel  atoms  with  4  layers.  There  are  27  atoms  in 
the  first  and  fourth  layers,  and  33  atoms  in  the  second  and  third  layers.  Figure  lb  illustrates  the  model 
used  in  this  calculation,  and  the  reactant,  product,  and  transition  state  sites  are  labelled  R,  P,  and  T, 
respectively.  Four  surface  atoms  (shaded  in  the  figure)  and  three  second-layer  atoms  surrounding  the 
reactant  site  were  allowed  to  move. 

Two  different  models  were  used  in  the  Ni(l  10)  calculations,  because  there  are  two  different  diffusional 
directions;  one  along  the  [llO]  direction  (parallel  to  the  channel)  and  the  other  along  the  [001]  direction 
(perpendicular  to  the  channel).  Hie  model  for  calculation  of  diffusion  in  the  [llO]  direction  consists  of 
186  atoms  over  3  layers  and  is  illustrated  in  Figure  1c.  There  are  42  atoms  in  the  first,  third,  and  fifth 
layers,  and  30  atoms  in  the  second  and  fourth  layers.  Six  surface  atoms  (darkly  shaded)  and  four  second- 
layer  atoms  (lightly  shaded)  adjacent  to  the  labelled  reactant,  product,  and  transition  state  sites  were 
allowed  to  move.  The  model  for  calculation  of  diffusion  in  the  [001]  direction  (Figure  Id,)  has 
219  atoms,  with  49  atoms  in  each  of  the  first,  third,  and  fifth  layers;  and  36  atoms  in  each  of  the  second 
and  fourth  layers.  Atoms  (six  in  the  first  layer  and  four  in  the  second  layer)  adjacent  to  reactant,  product, 
and  transition  state  sites  (labelled  R,  P,  and  T,  respectively)  were  allowed  to  move  and  are  shaded. 
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Figure  la.  Ni(lOO).  Figure  lb.  Ni(lll). 

Crystal  Lattice  Models  Used  in  the  Calculations  of  a  Nickel  Adsorbate 


Figure  lc.  Crystal  Lattice  Models  Used  in  the  Calculations  of  a  Nickel  Adsorbate  on  Ni(l  10)  (for 
calculations  of  diffusion  in  the  fllOl  direction) 


The  number  of  atoms  allowed  to  move  in  any  of  the  crystal  models  studied  was  not  varied  to  check 
for  convergence  of  a  diffusion  rate.  The  substrate  models  were  arbitrarily  chosen  and  used  in  each 
calculation  for  comparative  purposes  only.  However,  identical  lattice  models  were  used  for  all  the 
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Figure  Id.  Crystal  Lattice  Models  Used  in  the  Calculations  of  a  Nickel  Adsorbate  on  Ni(llO)  (for 
calculations  of  diffusion  in  the  TOOll  direction). 

potential  models  because  our  interest  was  not  in  calculating  the  absolute  true  rate,  because  of  the  lack  of 
consistent  experimental  data  with  which  to  compare,  but  rather  to  investigate  how  the  potentials  compared 
to  each  other.  Indeed,  for  estimating  a  very  precise  diffusion  rate,  it  may  be  necessary  to  allow  for  more 
substrate  atoms  to  relax  because  of  the  coupling  between  the  diffusing  Ni  atom  and  the  Ni -substrate  atom 
motion.  However,  in  conjunction  with  the  results  available  in  the  research  literature,  which  is  discussed 
in  Section  3.3.1,  the  calculated  rates  may  not  vary  significantly  when  more  substrate  atoms  are  allowed 
to  move. 

3.  RESULTS  AND  DISCUSSION 

Diffusion  coefficients  for  self-diffusion  of  nickel  on  the  (100)  and  (1 10)  surfaces  are  given  in  Table  2; 
Arrhenius  parameters  are  extracted  from  linear  least-squares  fits  of  In  [D(T)]  over  the  temperature  ranges 
listed  in  Table  2  and  are  given  in  Table  3,  along  with  reported  experimental  values.  The  tabulated 
diffusion  coefficients  as  a  function  of  temperature  are  provided  to  reflect  on  the  accuracy  of  fit  to 
Arrhenius  form  over  a  temperature  range  when  required. 

3.1  Surface  Structure  Effects.  Comparison  of  the  classical  diffusional  barriers  for  Ni/Ni(100), 
Ni/Ni(l  10),  and  Ni/Ni(l  1 1)  indicates  that  diffusion  will  occur  most  readily  on  die  (1 1 1)  surface,  and  more 
readily  along  the  channels  of  atoms  on  the  Ni(l  10)  surface,  which  is  typical  of  the  FCC  metals  (Ehrlich 


Table  2.  Calculated  Diffusion  Coefficients  (cm2/s) 


T(K) 

FBD* 

VC* 

or 

ATVF4 

Ni/Ni(100) 

50.0 

1.286  (-67) 

6.454  (-74) 

2.101  (-48) 

2.286  (-96) 

100.0 

1.786  (-35) 

1.270  (-38) 

7.185  (-26) 

8.863  (-50) 

200.0 

2.463  (-19) 

6.706  (-21) 

1.545  (-14) 

2.254  (-26) 

223.0 

1.145  (-17) 

4.548  (-19) 

2.295  (-13) 

5.931  (-24) 

243.0 

1.789  (-16) 

9.308  (-18) 

1.584  (-12) 

3.207  (-22) 

263.0 

1.840  (-15) 

1.204  (-16) 

8.154  (-12) 

9.464  (-21) 

283.0 

1.362  (-14) 

1.086  (-15) 

3.331  (-11) 

1.732  (-19) 

300.0 

6.054  (-14) 

5.589  (-15) 

9.510  (-11) 

1.513  (-18) 

400.0 

3.014  (-11) 

5.129  (-12) 

7.500  (-09) 

1.255  (-14) 

600.0 

1.504  (-08) 

4.724  (-09) 

5.936  (-07) 

1000.0 

2.170  (-06) 

1.113  (-06) 

1.965  (-05) 

1.452  (-07) 

1100.0 

4.276  (-06) 

2.346  (-06) 

3.167  (-05) 

3.893  (-07) 

1200.0 

7.524  (-06) 

4.364  (-06) 

4.714  (-05) 

8.860  (-07) 

1300.0 

1.214  (-05) 

7.381  (-06) 

6.601  (-05) 

1.777  (-06) 

1400.0 

1.829  (-05) 

1.158  (-05) 

8.809  (-05) 

3.226  (-06) 

1500.0 

1.711  (-05) 

1.131  (-04) 

5.410  (-06) 

1600.0 

3.560  (-05) 

2.407  (-05) 

1.408  (-04) 

8.504  (-06) 

Ni/Ni(U0)[ll0] 

50.0 

7.276  (-54) 

1.711  (-56) 

2.592  (-47) 

3.587  (-68) 

60.0 

1.628  (-45) 

1.062  (-47) 

4.720  (-40) 

1.971  (-57) 

80.0 

4.712  (-35) 

1.107  (-36) 

5.953  (-31) 

5.580  (-44) 

100.0 

9.263  (-29) 

4.714  (-30) 

1.788  (-25) 

6.871  (-36) 

125.0 

1.027  (-23) 

9.723  (-25) 

4.417  (-21) 

2.101  (-29) 

150.0 

2.394  (-20) 

3.435  (-21) 

3.791  (-18) 

4.500  (-25) 

200.0 

3.917  (-16) 

9.465  (-17) 

1.783  (-14) 

1.190  (-19) 

Note:  Power  of  10  m  parenthesis. 

*Foiles,  Btskes.  and  Daw  approach 
*Voter  and  Qten  approach 
cOh  and  Johnson  approach 
‘Ackiand,  Tichy,  Vitek,  and  Frnnis  approach 
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Table  2.  Calculated  Diffusion  Coefficients  (cm2/s)  (continued) 


T(K) 

FBD* 

VCb 

Of 

ATVF* 

250.0 

1.327  (-13) 

4.387  (-14) 

2.865  (-12) 

2.159  (-16) 

300.0 

6.459  (-12) 

2.631  (-12) 

8.489  (-11) 

3.227  (-14) 

400.0 

8.308  (-10) 

4.393  (-10) 

5.880  (-09) 

1.696  (-11) 

500.0 

1.532  (-08) 

9.471  (-09) 

7.485  (-08) 

7.299  (-10) 

600.0 

1.069  (-07) 

7.337  (-08) 

4.082  (-07) 

8.975  (-09) 

1,100.0 

8.844  (-06) 

7.689  (-06) 

1.930  (-05) 

2.701  (-06) 

1 ,200.0 

1.375  (-05) 

1.213  (-05) 

2.838  (-05) 

4.780  (-06) 

1,300.0 

1.998  (-05) 

1.799  (-05) 

3.934  (-05) 

7.750  (-06) 

1,400.0 

2.753  (-05) 

2.521  (-05) 

5.204  (-05) 

1.173  (-05) 

1,500.0 

3.633  (-05) 

3.377  (-05) 

6.633  (-05) 

1.679  (-05) 

1,600.0 

4.632  (-05) 

4.362  (-05) 

8.201  (-05) 

2.299  (-05) 

Ni/Ni(110)  [001] 

50.0 

2.992  (-134) 

3.346  (-171) 

1.387  (-104) 

8.352  (-157) 

60.0 

2.122  (-112) 

3.420  (-143) 

1.141  (-87) 

3.460  (-131) 

80.0 

4.611  (-85) 

3.731  (-108) 

1.689  (-66) 

3.876  (-99) 

100.0 

1.209  (-68) 

4.103  (-87) 

8.891  (-54) 

6.897  (-80) 

125.0 

1.690  (-55) 

2.875  (-70) 

1.378  (-43) 

1.794  (-64) 

150.0 

9.932  (-47) 

4.964  (-59) 

8.715  (-37) 

3.461  (-54) 

200.0 

9.228  (-36) 

5.629  (-45) 

2.826  (-28) 

2.554  (-41) 

250.0 

3.543  (-29) 

1.541  (-36) 

3.658  (-23) 

1.364  (-33) 

300.0 

8.716  (-25) 

6.522  (-31) 

9.412  (-20) 

1.947  (-28) 

400.0 

2.695  (-19) 

7.078  (-24) 

1.736  (-15) 

5.100  (-22) 

500.0 

5.166  (-16) 

1.181  (-19) 

3.215  (-13) 

2.913  (-18) 

600.0 

6.531  (-14) 

7.723  (-17) 

1.426  (-11) 

8.932  (-16) 

Note:  Power  of  10  in  parentheses. 

*Foiles,  Basko,  and  Daw  approach 
^Voter  and  Chen  approach 
cOh  md  Johnson  spproach 
dAckland,  Tichy,  Vitek  and  Fmnis  approach 


Table  3.  Arrhenius  Parameters  and  Static  Properties 


Theoretical 


FBD*  VCb 


Experimental 

FIM 

Mass  Transfer 

Tung  & 

Bonzel  &  Maiya  & 

Graham  Liu  et  al.e 

Latta  Blakely 

(1980)  (1991) 

(1978)  (1965) 

Ni/Ni(100) 


E,  (eV)  0.64 


D0  (cm2/s)  0.003 


0.64 


Az  (A)*  1.87 


0.70 

0.45 

0.93 

0.0037 

0.0035 

0.0064 

0.71 

0.45 

0.94 

1.83 

2.00 

1.93 

Ni/Ni(U0)[ll0] 


E.(eV) 

0.50 

0.53 

0.44 

0.64 

D0  (cm2/s) 

0.0033 

0.0038 

0.0036 

0.0045 

Vf  (eV)f 

0.50 

0.53 

0.44 

0.65 

Az  (A)* 

1.11 

1.06 

1.20 

1.14 

0.23 

0.45 

0.76 

1.85 

10'9 

0.001 

0.009 

23.9 

Ni/Ni(110)[001] 


E.(eV) 

1.30 

1.67 

1.01 

1.53 

D0  (cm2/s) 

0.0065 

0.0075 

0.0065 

0.0077 

Vf  (eV)f 

1.31 

1.68 

1.02 

1.54 

Az  (A)* 

1.11 

1.06 

1.20 

1.14 

>.32 

0.45 

1.95 

1.74 

10'7 

0.001 

470.0 

12.8 

Ni/Ni(lll) 


E,(eV) 


Az  (A)* 


0.06 

0.07 

0.05 

0.05 

1.84 

1.79 

1.94 

1.89 

•Foiles,  Basket,  and  Daw  approach 

^Voter  and  Chen  approach 

cOh  and  Johnson  approach 

dAchland,  Tichy,  Vitek,  and  Fiimis  approach 

•Reanalyzed  data  of  Tung  and  Graham  (1980)  (see  text) 

‘Classical  diffusicmal  barrier 

*  Distance  above  the  surface  in  die  reactant  binding  site 


13 


and  Stolt  1980;  Doll  and  Voter  1987;  Bonzel  1983;  Ehrlich  1990;  De  Lorenzi  1986;  De  Lorenzi  and 
Jacucci  1985).  Additionally,  all  of  the  models  indicate  two  distinct  binding  sites  on  the  (111)  surface, 
which  was  observed  in  experiment  (Tung  and  Graham  1980).  Recent  lattice  statics  calculations  (Liu, 
Cohen,  Adams,  and  Voter  1991)  for  all  the  EAM-type  FCC  metals  have  led  to  similar  results. 

The  variations  in  diffusion  rates  have  been  explained  in  terms  of  the  influence  of  atomic  packing  in 
the  different  surface  structures.  The  FCC  (100)  surface  is  loosely  packed  compared  to  the  more  compact 
(111)  surface.  Similarly,  in  the  [1 10]  direction,  the  atomic  structure  of  the  (1 10)  surface  is  more  compact 
relative  to  that  in  die  [001]  direction.  Detailed  MD  studies  (De  Lorenzi,  Jacucci,  and  Pontikis  1982), 
using  Lennard-Jones  interatomic  potentials  for  surfaces  of  rare  gas  solids,  have  shown  clearly  that  close 
packing  enhances  adatom  mobility.  Our  results,  together  with  the  recent  calculations  (Liu,  Cohen,  Adams, 
and  Voter  1991)  for  the  FCC  metals  described  by  EAM  models,  confirm  these  findings.  Further  details 
of  surface  structure  effects  and  the  pertinent  adatom  diffusion  mechanisms  are  discussed  in  Section  3.3. 

3.2  Embedded-Atom  Functional  Effects.  Comparison  of  diffusion  rates  among  the  four  EAM 
potential  energy  models  shows  that  in  all  cases,  the  OJ  model  predicts  the  fastest  diffusion  rates,  followed 
by  the  FBD  model.  The  VC  model  predicts  a  faster  diffusion  rate  than  the  ATVF  model  for  diffusion  on 
the  (100)  surface  and  in  the  [110]  direction  on  the  (110)  surface,  but  the  ordering  of  diffusion  rates  is 
switched  for  diffusion  in  the  [001]  direction  on  the  (1 10)  surface.  Dependence  of  atom  diffusion  rates  on 
the  pair  functionals  has  not  been  analyzed  to  date.  In  order  to  analyze  the  different  diffusion  rates  based 
on  different  EAM  models,  we  first  examined  the  component  terms;  such  as  p(r),  F(p),  and  <D(r),  that 
constitute  an  embedded-atom  pair  functional.  Furthermore,  such  an  examination  is  useful  in  its  own  right 
because  of  the  different  functional  forms  and  of  the  empirical  procedures  employed  in  devising  the  various 
EAM  models  as  discussed  briefly  in  Section  2.2. 

The  electron  density  variation,  p(r),  in  the  modeled  nickel  is  characterized  differently  in  each  of  the 
adapted  pair  functionals.  The  electron  density  is  a  simple  superposition  of  contributions  from  neighboring 
atoms.  In  order  to  facilitate  a  suitable  comparison  of  p(r)  and  the  resulting  F(p)  of  the  different  EAM 
models,  p(r)  is  scaled  to  its  equilibrium  value,  pe,  the  electron  density  of  an  atom  in  the  bulk  metal  at 
equilibrium.  Such  a  normalized  p(r)  for  all  the  used  EAM  models  is  shown  in  Figure  2.  Excepting  for 
the  ATVF  model,  all  the  other  functionals  show  exponential  behavior.  The  equilibrium  electron  densities 
are  0.067,  0.38,  12.55,  and  37.41  A'3  for  FBD,  VC,  OJ,  and  ATVF  models,  respectively.  The  targe 
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FBD 

VC 


Figure  2.  Normalized  Electron  Density  Plots  for  Nickel  as  a  Function  of  Interatomic  Distance:  oc  is 
the  Atomic  Electron  Density  in  the  Bulk  Metal  at  Equilibrium. 
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differences  in  pe  and  p(r)  are  due  to  the  different  procedures  used  for  mimicking  the  solid-state  bulk  metal 
atom  effects. 

The  embedding  functions  based  on  the  respective  p(r)’s  are  shown  in  Figure  3.  Excepting  for  the  VC 
model,  for  which  the  pair  potential  is  taken  to  be  of  Morse  type,  the  other  embedding  function  shapes  are 
similar  and  the  functions  are  attractive  in  nature  as  the  corresponding  pair  interactions  are  treated  to  be 
repulsive.  Comparison  of  diffusion  rates  based  on  VC  model  with  that  obtained  by  using  other  EAM 
models  indicate  that  the  two-body  potential  can  be  attractive  or  repulsive,  with  the  corresponding 
embedding  function  being  repulsive  or  attractive,  respectively,  without  much  effect  on  calculated  diffusion 
rates.  This  may  be  reasoned  out  as  due  to  the  following:  the  recipe  of  fitting  to  the  equation  of  state 
brings  about  some  kind  of  balance  between  the  two  terms.  Furthermore,  the  function  F(p)  and  d>(r)  need 
not  be  uniquely  determined,  as  discussed  by  Daw  and  Baskes  (1983)  on  EAM  approach,  by  the  empirical 
procedures.  As  the  examination  of  the  component  terms,  F(p)  and  <h(r),  did  not  lead  to  a  physically 
plausible  explanation  for  the  variation  in  diffusion  rates  based  on  the  different  EAM  models,  the  next  step 
we  chose  was  to  examine  the  effective  two-  and  three-body  potentials  that  can  be  constructed  from  an 
EAM  model. 

Based  on  detailed  MD  studies  of  surface  diffusion  on  FCC  and  body-centered  cubic  (BCC)  systems, 
Jacucci  and  coworkers  have  shown  that  the  rate  of  diffusion  is  expected  to  be  affected  by  the  degree  of 
stiffness  in  the  repulsive  wall  (De  Lorenzi  and  Jacucci  1985).  To  investigate  this,  following  Foiles  (1985, 
1987),  we  approximated  the  EAM  energy  as  an  effective  pair  potential  by  expanding  the  embedding 
functional  in  a  Taylor  series  and  dropping  three-body  and  higher  order  terms.  We  have  thus  calculated 
the  effective  pair  potential  OeR(r)  as 

*eff (0  *  ®<r)  ♦  2F'  (pe)  p(r)  ♦  F"  (pe)  [  p(r)  ]2  (18) 

as  a  function  of  interatomic  distance  r  for  each  of  the  four  EAM  potential  functions  described  above  and 
show  the  effective  potentials  in  Figure  4.  F(Pj)  and  F'ipJ  in  Equation  18  are  the  first  and  second 
derivatives  of  the  embedding  energy  with  respect  to  the  local  election  density,  p. 

The  effective  pair  potential  curve  for  the  ATVF  model  differs  from  die  other  curves  for  the  FBD,  VC, 
and  OJ  models  both  in  the  repulsive  wall  region  (located  at  a  smaller  R/a  value)  and  at  R/a  values 
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greater  than  the  minimum  of  the  well  where  a  "bump"  exists  rather  than  merely  an  attractive  tail. 
However,  on  the  basis  of  simulations  on  self-  and  defect-diffusion  studies  in  bulk  and  in  grain  boundaries, 
one  may  safely  conclude  that  the  tail  region  effects  are  small  for  adatom  diffusion.  Although  the  ATVF 
potential  energy  curve  has  a  well  depth  approximately  equal  to  the  OJ  curve,  its  repulsive  wall  is  much 
softer  than  the  other  effective  pair  potentials.  From  Table  3  and  Figure  5,  we  find  that  the  ATVF  model 
consistently  predicts  larger  barriers  for  diffusion  (correspondingly  lower  diffusivities)  on  the  (100)  surface 
and  in  the  [110]  direction  on  the  (1 10)  surface.  The  predicted  large  diffusion  barriers  predicted  by  the 
ATVF  model  thus  concur  with  the  conclusion  from  the  detailed  MD  studies  (De  Lorenzi  and  Jacucci  1985) 
that  softness  hinders  adatom  mobility. 

The  effective  potential  energy  curves  for  the  FBD,  VC,  and  OJ  models  have  a  similar  shape;  the  main 
differences  occur  in  the  potential  well  depth  and  in  the  long-range  attractive  tails.  While  OJ  and  FBD 
curves  are  superimposable  up  to  the  well  depth,  compared  to  OJ  and  FBD  curves,  the  VC  curve  differs 
from  the  repulsive  region  itself.  Of  the  three  effective  potentials  which  have  the  same  shape,  the  OJ 
potential  energy  curve  has  the  deepest  attractive  well,  and  the  FBD  curve  has  the  next  deepest  well.  The 
activation  energies  in  Table  3  and  the  Arrhenius  plots  in  Figure  5  indicate  that  the  diffusion  rates  predicted 
by  these  models  are  also  in  the  same  order  in  terms  of  fastest  to  slowest  and  low  to  high  with  respect  to 
barrier  height 

Effect  of  higher-order  terms  in  the  Taylor  series  of  F(p)  on  d>efJ(r)  is  evident  from  a  comparison  with 
the  inset  in  Figure  2,  which  corresponds  to  the  <I>en(r)  computed  from  only  the  first  two  terms  in 
Equation  18.  Furthermore,  the  inclusion  of  three-body  interactions,  which  is  more  cumbersome  to 
evaluate,  and  higher-order  terms  may  lead  to  a  more  definitive  analysis,  in  terms  of  the  effective 
potentials,  for  the  diffusion  rates  predicted  by  FBD,  VC,  and  OJ  models. 

3.3  Comparison  With  Experiment.  Experimental  results  presented  are  from  FIM  and  mass  transfer 
experiments;  however,  the  results  shown  in  Table  3  indicate  that  contradictory  experimental  data  for  self¬ 
diffusion  on  nickel  exist.  Therefore,  stringent  comparison  of  the  predicted  result  to  experiment  is  difficult 
due  to  the  wide  variation  of  experimental  results  shown  in  Table  3.  FIM  measurements  provide  rates  of 
intrinsic  diffusion,  and  are  essentially  independent  of  surface  defects,  while  mass  transfer  measurements 
are  performed  over  larger  distances  on  die  crystal,  thus  involving  a  range  of  orientations  and  surface 
defects.  The  rates  of  diffusion  in  our  calculations  would,  therefore,  most  likely  correspond  to  the  intrinsic 
diffusion  rates  measured  by  FIM. 
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Effective  Pair  Potential  Energy  Curves  as  a  Function  of  Interatomic  Distance  Scaled  to  the 
Lattice  Parameter  of  Nickel.  3.52  X. 
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Figure  5.  Arrhenius  Plot  of  Self-Diffusion  on  (100)  Surface  and  on  (1 10)  Surface  in  the  Directions  of 
[lfOl  and  fOOll:  (solid  line)  FBD,  (long-dashed  line)  VC.  (short-dashed  line)  OJ.  and 
(dot-dashed  line)  ATVF  model  predictions. 
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3.3.1  The  (100)  Surface.  Tung  and  Graham  (1980),  in  a  discussion  of  their  FIM  results,  indicate  that 
while  reliable  diffusion  measurements  for  the  two-dimensional  (100)  and  (111)  surfaces  are  difficult  to 
obtain,  the  onset  temperature  for  adatom  mobility  on  these  surfaces  can  be  accurately  determined.  By 
determining  the  onset  temperature,  assigning  a  value  of  D  =  1017  cm2/s  at  the  onset  temperature,  and 
assuming  a  preexponential  factor  of 


kBTl 
2  h 


2 


(19) 


where  1  denotes  tire  hop  length  between  binding  sites,  Tung  and  Graham  (1980)  estimate  activation 
energies  of  0.63  and  0.33  eV  for  self-diffusion  on  the  (100)  and  (111)  surfaces,  respectively.  If  the 
assumptions  made  by  Tung  and  Graham  (1980)  are  reasonable,  all  models  are  in  good  agreement  with 
experiment  for  diffusion  on  the  (100)  surface,  but  agreement  is  poor  for  the  Ni(lll)  surface.  Recent 
lattice  statics  calculations  (Liu,  Cohen,  Adams,  and  Voter  1991)  using  the  VC  and  another  set  (Adams, 
Foiles,  and  Wolfer  1989)  of  FBD-type  models  have  led  to  values  of  0.68  and  0.63  eV  for  the  activation 
energy  of  self-diffusion  on  (100)  surface.  These  are  in  good  agreement  with  our  estimated  values  of  0.71 
and  0.64  eV  for  the  VC  and  FBD  models,  respectively.  The  respective  preexponential  factors,  Dq,  5.4, 
and  1.6  1(T3  cm2/s  are  comparable  to  our  values,  3.7  and  3.5  10'3  cm2/s.  These  results  are  based  on  the 
hopping  mechanism.  Discussion  involving  energetics  for  more  complicated  mechanisms  such  as  adatom 
exchange  (Bassett  and  Webber  1978;  Halicioglu  and  Pound  1979)  is  included  in  Section  3.3.3. 


The  agreement  for  the  energetics  based  on  hopping  mechanism  between  Liu,  Cohen,  Adams,  and 
Voter’s  (1991)  approach  and  our  VTST  calculations  (which  allowed  only  a  small  number  of  the  substrate 
atoms  w  move)  is  remarkable  and  indicates  that  the  convergence  of  our  estimated  diffusion  rates  is  more 
than  adequate  for  the  present  comparative  studies  of  potential  model  effects.  The  results  of  these 
potential-based  studies  are  also  in  accord  with  the  recent  first-principles  calculations  (Feibelman  1990)  of 
self-diffusion  of  A1  on  Al(001).  The  latter  studies  have  discussed  that  the  atomic  relaxations  in  A1  have 
a  small  effect  on  the  adatom  diffusion  barrier  height.  Liu,  Cohen,  Adams  and  Voter  (1991)  found  that 
typically  two  to  six  moving  atom  layers  resulted  in  activation  energy  convergence  to  less  than  0.01  eV. 
These  results  in  conjunction  with  ours  would  imply  that  self-diffusion  on  the  relaxed  surface  would  be 
slightly  different  from  a  rigid  substrate  and  that  effects  of  surface  motion  on  rates  would  be  minimal. 


3.3.2  The  (111)  Surface.  The  reliability  of  the  estimated  activation  energy  of  0.33  eV  for  self¬ 
diffusion  on  the  (111)  surface  is  questionable  due  to  the  existence  of  two  distinct  binding  sites  on  this 
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surface,  and  Ae  inability  of  Tung  and  Graham  (1980)  to  determine  the  onset  temperature  for  adatom 
displacement  between  the  two  sites.  The  onset  temperature  used  for  the  estimate  of  the  activation  energy 
for  diffusion  on  the  (111)  surface  corresponds  to  movements  of  adatoms  over  approximately  two  lattice 
spacings,  and  does  not  correspond  to  motion  from  the  two  binding  sites  shown  in  Figure  lb.  All  of  the 
four  EAM  potentials  predict  two  distinct  binding  sites  (denoted  in  Figure  lb  as  R  and  P)  in  good 
agreement  with  experiment  (Tung  and  Graham  1980),  but  the  differences  in  energies  for  the  two  binding 
sites  are  less  than  0.03  eV  for  all  four  potentials,  and  the  barriers  to  diffusion  between  the  two  sites  are 
no  greater  than  0.07  eV  for  all  four  EAM  potentials  (see  Table  3).  Rates  of  self -diffusion  on  Ni(l  1 1)  were 
not  calculated  because  when  barriers  are  too  low,  die  approximation  of  uncorrelated  hops  assumed  in 
Equation  1  becomes  invalid  and  other  approaches  for  calculating  diffusion  coefficients  are  needed. 

3.3.3  The  (110)  Surface  and  Diffusion  Anisotropy.  Comparison  of  the  calculated  diffusion  rates  with 
experimental  values  on  the  Ni(110)  surface  is  also  not  straightforward  because  the  experimental  results 
are  questionable  (Bonzel  1983).  The  data  provided  by  Tung  and  Graham  (1980)  for  diffusion  in  the  [110] 
direction  indicate  an  activation  energy  of  0.23  eV  and  an  extremely  small  preexponential  factor.  Rates 
of  diffusion  measured  from  FIM  experiments  typically  tend  to  be  much  larger  than  those  measured  from 
mass  transfer  methods  (by  at  least  10*)  (Bonzel  1983).  The  FIM  data  for  Ni(l  10)  presented  by  Tung  and 
Graham  (1980),  however,  appear  to  be  at  least  104  smaller  than  die  corresponding  mass  transfer 
measurement,  in  contradiction  to  typical  comparisons  of  FIM  (Tung  and  Graham  1980)  results  to  mass 
transfer  results  (Bonzel  and  Latta  1978;  Maiya  and  Blakely  1965).  This  has  led  Bonzel  (1983)  in  his 
review  of  nickel  self-diffusion  results  to  question  the  reliability  of  the  FIM  measurements  for  clean  nickel. 
Reccndy,  Liu,  Cohen,  Adams,  and  Voter  (1991),  in  connection  with  their  calculations  of  surface  diffusion 
based  on  EAM  models,  have  expressed  some  concern  about  fitting  an  Arrhenius  form  to  a  few  FIM  data 
points  and  the  resulting  large  errors  in  estimates  of  Dq.  Both  our  estimates,  as  well  as  their  estimates,  are 
about  1(T3  cm2/s  compared  to  the  value  of  10'9  cm2/s  reported  by  Tung  and  Graham.  Liu,  Cohen,  Adams, 
and  Voter  (1991)  have  thus  re-fit  the  Tung  and  Graham  (1980)  data  to  an  Arrhenius  equation,  assuming 
a  preexponential  factor  of  0.001  cm2/s.  The  reanalyzed  Arrhenius  parameters  are  shown  in  Table  3.  The 
OJ,  FBD,  and  VC  models  are  in  reasonable  agreement  with  the  reanalyzed  experimental  result  for 
diffusion  along  the  [1 10]  direction.  The  reanalyzed  diffusion  rates,  however,  present  discrepancies  for 
diffusion  along  die  [001]  direction. 

The  first  analysis  of  the  Tung  and  Graham  (1980)  data  (using  D0  =  10'9  cm2/s)  suggests  directional 
anisotropy  in  diffusion  on  the  Ni(U0)  surface;  reanalysis  of  the  data  (D0  =  10‘3  cm2/s)  removes  the 
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directional  anisotropy  of  diffusion  on  the  Ni(l  10)  surface  and  suggest  an  activation  energy  of  0.45  eV  for 
diffusion  along  both  directions.  Directional  anisotropy  on  Ni(l  10),  however,  is  clearly  evident  in  the  mass 
transfer  results  of  Bonzel  and  Latta  (1978).  To  confuse  the  issue  further,  directional  anisotropy  on  Ni(l  10) 
is  not  evident  in  mass  transfer  studies  by  Maiya  and  Blakely  (1965);  in  fact,  no  effect  of  surface 
orientation  on  diffusion  rates  is  observed  for  the  (100),  (111),  and  (110)  surfaces,  leading  Bonzel  (1983) 
in  his  review  to  suggest  that  the  nickel  surfaces  used  in  the  study  were  insufficiently  clean. 

All  four  EAM  models  predict  a  directional  anisotropy,  and  that  adatom  motion  parallel  to  the  channels 
is  more  probable  than  hopping  between  channels.  This  is  clearly  borne  out  in  the  Arrhenius  plots 
(Figure  5)  for  the  diffusion  on  (110)  surface  in  [llO]  and  [001]  directions.  Our  calculated  diffusion 
coefficients  in  the  [001]  direction  are  closer  to  the  mass  transfer  results  than  the  FIM  data  given  by  Tung 
and  Graham  (1980)  or  reanalyzed  by  Liu,  Cohen,  Adams,  and  Voter  (1991);  however,  Bonzel  (1983) 
suggests  that  because  of  the  large  preexponential  factor  attained  from  the  mass  transfer  experiments,  these 
Arrhenius  parameters  correspond  to  a  nonlocalized  diffusion  process  (which  is  by  its  nature  isotropic),  and, 
thus,  the  true  diffusion  rate  in  the  [001]  direction  is  not  measured. 

Diffusion  across  the  channels  on  the  (1 10)  surface  via  adatom-exchange  with  one  of  the  channel  atoms 
has  been  observed  (Wrigley  and  Ehrlich  1980).  Additionally,  there  is  evidence  from  the  recent  first- 
principles  total-energy  calculations  that  this  complex  mechanism  for  diffusion  can  be  invoked  for  surface 
self-diffusion  on  the  (100)  surface  (Kellog  and  Feibelman  1990).  Liu,  Cohen,  Adams,  and  Voter  (1991) 
calculate  Arrhenius  parameters  for  this  atom  exchange  mechanism  for  diffusion  across  the  channels  on  the 
Ni(110)  surface,  which  are  in  reasonable  agreement  with  the  reanalyzed  Tung  and  Graham  (1980)  data, 
and  which  predict  a  much  lower  activation  energy  (0.39  and  0.42  eV  for  VC  and  FBD-type  models)  for 
diffusion  via  adatom-exchange  compared  to  our  values  (0.53  and  0.50  eV  for  VC  and  FBD  models)  based 
on  hopping  across  the  channels.  The  respective  D0’s,  4.0  and  2.3  x  10'3  cm2/s,  are  in  agreement  with  our 
values,  3.8  and  3.3  x  10'3  cm2/s. 

Recently,  Hansen,  Stoltze,  Jacobsen,  and  N0rskov  (1991)  examined  diffusion  on  the  (100)  and  (110) 
surfaces  of  copper  via  adatom-exchange  and  found  that  the  energies  calculated  from  the  effective  medium 
theory  predict  the  adatom  exchange  mechanism  as  the  energetically  most  favorable  mechanism  for 
diffusion  on  these  surfaces.  Furthermore,  they  found  that  a  simple  hopping  mechanism  is  energetically 
favored  for  self-diffusion  on  (1 1 1)  surface  and  in  the  [001]  direction  on  the  (1 10)  surface.  While  these 
results  thus  indicate  directional  anisotropy  in  diffusion  on  the  Cu(l  10)  surface,  there  appears  to  be  some 
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disagreement  between  their  findings  and  Liu,  Cohen,  Adams,  and  Voter  (1991)  findings  concerning  the 
most  energetically  favorable  mechanism  for  the  Cu(100)  surface  as  well.  The  latter  studies  based  on  VC 
and  FB EX- type  EAM  models  show  that  Ni  and  Cu  do  not  favor  adatom  exchange  mechanism  for  diffusion 
on  (100)  surface  and  that  self-diffusion  on  (110)  surface,  based  on  the  exchange  mechanism,  is  isotropic. 
Although  we  did  not  calculate  rates  of  diffusion  invoking  this  mechanism  on  either  the  (100)  and  (1 10) 
surfaces,  our  results,  when  compared  to  the  Liu,  Cohen,  Adams,  and  Voter  (1991)  results,  would  suggest 
that  the  adatom  exchange  mechanism  describes  diffusion  in  the  [001]  direction  on  the  (110)  surface. 
Detailed  MD  studies  (De  Lorenzi,  Jacucci,  and  Pontikis  1982)  have  shown  that  the  predominance  of  the 
adatom  exchange  mechanism  results  in  nearly  isotropic  diffusion  which  is  further  confirmed  by  recent 
calculations  of  Liu,  Cbhen,  Adams,  and  Voter  (1991)  and  that  the  superposition  of  all  the  other  possible 
mechanisms  leads  to  diffusion  anisotropy  on  the  (110)  surface. 

4.  SUMMARY  AND  CONCLUSION 

Calculations  of  Ni  surface  self-diffusion  based  upon  VTST  were  performed  using  different  EAM 
potentials  in  an  effort  to  analyze  their  suitability  in  calculating  surface  and  interfacial  properties.  In  spite 
of  the  difficulty  of  comparing  the  calculated  results  with  the  measured  values,  the  EAM  potentials  were 
successful  in  reproducing  several  observables  of  nickel  adsorbed  on  its  surfaces.  On  the  basis  of  the 
hopping  mechanism,  all  of  the  models  predict  the  directional  anisotropy  for  self-diffusion  on  the  Ni(l  10) 
surface,  that  diffusion  occurs  more  easily  in  the  [110]  direction,  and  that  two  distinct  binding  sites  of  Ni 
exist  on  the  (1 1 1)  surface.  Estimated  classical  diffusion  barriers  indicate  that  diffusion  occurs  most  easily 
on  the  (1 1 1)  surface  compared  to  other  surfaces.  At  high  temperatures,  the  models  agree  with  each  other, 
but  as  temperatures  decrease,  the  predicted  diffusivities  fall  within  a  scatterband  as  is  clear  from  Figure  5. 

Although  all  four  embedded-atom  pair  functionals  provide  Arrhenius  parameters  which  are  in 
reasonable  agreement  with  the  measured  experimental  values  for  &  If-diffusion  on  Ni(100),  the 
experimental  differences  in  the  rates  of  diffusion  on  the  other  surfaces  must  be  resolved  before  a  stringent 
comparison  with  the  rates  will  indicate  if  further  refinement  is  necessary  for  the  EAM  models.  It  is 
encouraging,  however,  that  each  of  these  models  predict  surface  properties  which  have  been  determined, 
and  predict  reasonable  behavior  of  diffusion  on  all  of  the  surfaces.  These  studies  involving  the  surface 
diffusivity  and  surface  structure  with  realistic  interaction  potentials  combined  with  recent  advances  in  first- 
principles  calculations  as  well  as  the  explosive  developments  in  the  thin-film  growth  techniques  should 
result  in  renewed  impetus  for  more  refined  measurements  of  adatom  diffusion. 
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